% izh models with no conductance integration for injecting conductances or
% currents (last param.)

function [spikes Gampa]=izhikevich(stim_conductance_exc, stim_conductance_inh,neuron,g_true_I_false)

global v

if ~exist('g_true_I_false')
    g_true_I_false = true;
end
if ~exist('neuron')
    neuron = 1;
end


switch neuron
    case 1
        vr=-60;
        vt=-40;
        k=3;
        C=100;
        a=0.01;
        b=5;
        c=-50;
        d=100;
        vp=35;
    case 2
        % pyramidal from Izh & Edelman, 2008
        vr=-60;
        vt=-50
        k=3;
        C=80
        a=0.01;
        b=5;
        c=-60;
        d=10;
        vp=50;
    case 3
        % fs cell.
        vr=-55;
        vt=-40;
        k=1;
        C=20;
        a=.15;
        b=8;
        c=-55;
        d=200;
        vp=25;
    case 4
        % TC cell from ibid., assuming only tonic firing model
        vr= -60;
        vt= -50;
        k=   1.6;
        C= 200;
        a=0.1;
        b=   0; % Assuming it is in tonic mode all the time.
        c=   -60;
        d=   10;
        vp=  40;
    case 5
        % chattering cell from Izhikevich book, Dynamical Systems in
        % Neuroscience
        vr=-60;
        vt=-40;
        k=1.5;
        C=50;
        a=0.03;
        b=1;
        c=-40;
        d=150;
        vp=25;

end




T=1000;

dt_ms=1;


spikes = 0;
n=round(T/dt_ms);
v=vr*ones(1,n); u=0*v;
si= stim_conductance_inh .*[ones(1,n)];
se= stim_conductance_exc .*[ones(1,n)];

% model conductances.
tau_ampa  = exp(-dt_ms/5);
tau_nmda  = exp(-dt_ms/150);
tau_gabaa = exp(-dt_ms/6);
tau_gabab = exp(-dt_ms/150);

E_AMPA	= 0.0;
E_NMDA	= 0.0;
E_GABAa	= -70;
E_GABAb = -90;
Ggabaa = 0;Ggabab = 0; Gampa = 0;Gnmda = 0;

Ggabaa=zeros(1,n);
Ggabab=zeros(1,n);
Gampa=zeros(1,n);
Gnmda=zeros(1,n);
I=zeros(1,n);

gain_ampa=1;
gain_nmda=1;
gain_gabaa=1;
gain_gabab=1;
% end model conductances.


for i = 1:n-1
    
    % conductances from synapses.
    %Ggabaa(i+1) = Ggabaa(i)*tau_gabaa + dt_ms*(1-tau_gabaa)*gain_gabaa*si(i);
    %Ggabab(i+1) = Ggabab(i)*tau_gabab + dt_ms*(1-tau_gabab)*gain_gabab*si(i);
    
    %Gampa(i+1) = Gampa(i)*tau_ampa + dt_ms*(1-tau_ampa)*gain_ampa*se(i);
    %Gnmda(i+1) = Gnmda(i)*tau_nmda + dt_ms*(1-tau_nmda)*gain_nmda*se(i);

    Ggabaa(i+1) = si(i);
    
    Gampa(i+1) = se(i);
    
    
    xx = (-80-v(i))/60;
	xx = xx.*xx;
	NMDAgate = xx./(1+xx);
    
     
    %g= ( Gampa(i+1)       +NMDAgate.*Gnmda(i+1)        + Ggabaa(i+1)         + Ggabab(i+1));
	%E= ( Gampa(i+1)*E_AMPA+NMDAgate.*Gnmda(i+1)*E_NMDA + Ggabaa(i+1)*E_GABAa + Ggabab(i+1)*E_GABAb);
    %I(i) = v(i).*g-E;
    if g_true_I_false
        I(i) = (v(i)-E_AMPA)*Gampa(i+1) + (v(i)-E_NMDA)*NMDAgate*Gnmda(i+1)+ ...
            (v(i)-E_GABAa)*Ggabaa(i+1)+(v(i)-E_GABAb)*Ggabab(i+1);
    else
        I(i) = se(i);
    end
    
    v1=v(i); u1=u(i);
    
    v1=v1+0.5*(k.*(v1-vr).*(v1-vt)-u1-I(i))./C;    % for numerical 
    %v(v>104)=104;   % Nothing higher than Eca2+
    v1(v1>51)=51;   % Nothing higher than typical action potential
    v1(v1<-90)=-90;  % cheating on numerical integration. Nothing lower than EK

    u1=u1+0.5*a.*(b.*(v1-vr)-u1);                   % step is 0.5 ms
 
    v1=v1+0.5*(k.*(v1-vr).*(v1-vt)-u1-I(i))./C;    % for numerical 
    %v(v>104)=104;   % Nothing higher than Eca2+
    v1(v1>51)=51;   % Nothing higher than Eca2+
    v1(v1<-90)=-90;  % cheating on numerical integration. Nothing lower than EK

    u1=u1+0.5*a.*(b.*(v1-vr)-u1);                   % step is 0.5 ms

    v(i+1)=v1; u(i+1)=u1;
    
    %v(i+1)=v(i)+.5*dt_ms*(k*(v(i)-vr)*(v(i)-vt)-u(i)-I(i))/C;
    %v(i+1)=v(i+1)+.5*dt_ms*(k*(v(i+1)-vr)*(v(i+1)-vt)-u(i)-I(i))/C;
    %u(i+1)=u(i)+dt_ms*a*(b*(v(i)-vr)-u(i));
    if v(i+1)>=vp
        v(i)=vp;
        v(i+1)=c;
        u(i+1)=u(i+1)+d;
        spikes = spikes + 1;
    end
end
figure(1)
subplot(3,1,1)
plot(dt_ms*(1:n),[v' u']);
subplot(3,1,2)
plot(dt_ms*(1:n),[Gampa' Gnmda' Ggabaa' Ggabab']);
title('conductances');
subplot(3,1,3);
plot(dt_ms*(1:n),I);
title('Isyn (pA)');
drawnow
